In den Geowissenschaften ist die Fernerkundung die einzige Messtechnik, die eine vollständige meßtechnische Abdeckung auch globaler Skalen ermöglicht. Zur erfolgreichen Nutzung gehören die souverände Anwendung von existierenden und notwendigerweise auch die Anpassung dieser bzw.die Entwicklung eigener Methoden.
Ein wichtiger Anwendungsfall in der Geo- oder Umweltinformatik ist die Erfassung von Veränderungen mittels Satelliten-, Flugzeug- und Drohnenbildern (Change Detectin Analysis). Oft wird dies in Verbindung mit bio- und geophysikalischen oder auch vom Menschen verursachten Prozessen genutzt, um ein tieferes Verständnis und die Möglichkeit zur Entwicklung von Vorhersagemodellen zu erhalten. Dabei sind Bildanalysemethoden von hervorgehobener Bedeutung, um Informationen zu extrahieren, die es erlauben, die zugrunde liegenden Prozesse zu identifizieren. Ein ständig wachsender Bereich ist die Integration von immensen Datenmengen sog. Big Data in die Analysen.
Die nahezu exponetiell wachsende und bereits gegenwärtig kaum noch zu bewältigende Flut an verfügbaren Bild- und Fernerkundungsdaten muss einfach und automatisiert zugänglich sein um sowohl für den wissenschaftlichen Erkenntnisgewinn als auch für gesellschaftliche Zukunftsaufgaben nutzbar zu sein. Wir beginnen mit einer anwendungs orientierten Aufgabe in dieser Übung.
Unverarbeitete Satellitenbilder sind nicht zwingend informativ. Unser Auge kann zwar ein Echtfarbenbild relativ schlüssig und intuitiv interpretieren, aber eine zuverlässige und reproduzierbare wissenschaftliche Interpretation erfordert andere Ansätze. Zudem können Bild analysemethoden zusätzliche,unsichtbare und spezifischere Informationen ableiten. Wir haben bereits einfache Indizes berechnet. Auch die Ableitung des NDVI oder der Oberflächenalbedo ist eine physikalisch begründete Umwandlung von Bildsignalen in eine Meßgröße.
Um nützliche oder zielführende Informationen, z. B. über die Bodenbedeckung in einem Gebiet, zu erhalten, müssen wir die Daten daher fragestellungszentriert analysieren. Der wohl bekannteste und gängigste Ansatz ist die überwachte Klassifizierung von Bilddaten in Kategorien die von Interesse sind.
Diese Übung führt Sie in die Klassifizierung von Satelliten- und Luftvermessungsdaten in R ein.
Wir werden uns mit den folgenden Themen beschäftigen:
Bitte beachten Sie, dass alle Arten der Klassifizierung eine ind der Regel aufwendige Datenvorverarbeitung erfordern. Im Zentrum stehen dann Modellbildung und Qualitätsschätzung die als handwerkliche Grundlagen der Klassifizierung betrachtet werden können um schließlich in dr Datennachverarbeitung als wesentlichen Bestandteil die inhaltliche Interpretation der Ergebnisse durchzuführen. Wir werden versuchen die notendigen Minimalanforderungen in der Kurseinheiten zu erarbeiten.
Bei der überwachten Klassifizierung von Landbedeckungen wird aus einer begrenzten Menge Trainingsdaten zur Landbedeckung ein Modell abgeleitet, das die jeweilige Landbedeckung für den gesamten Datensatz vorhersagt. Die Landbedeckungstypen werden also a priori definiert, und das Modell versucht, diese Typen auf der Grundlage der Ähnlichkeit zwischen den Eigenschaften der Trainingsdaten und dem Rest des Datensatzes vorherzusagen.

Ganz pragmatisch erfordern Klassifizierungsaufgaben im Allgemeinen die folgenden Schritte:
Die folgende Abbildung zeigt die Schritte einer überwachten Klassifikation im Detail. Die optionalen Segmentierungsoperationen sind obligatorisch für objektorientierte Klassifizierungen, die sich nicht nur auf die Geometrie der Objekte, sondern auch auf die Werte jeder einzelnen Rasterzelle im Eingabedatensatz stützen. Um einzelne Objekte wie Häuser oder Baumkronen abzugrenzen, verwenden Fernerkundungsexperten Segmentierungsalgorithmen, die die Homogenität der Pixelwerte innerhalb ihrer räumlichen Nachbarschaft berücksichtigen.

In diesem Tutorium werden die Sentinel-2-Bilder aus der vorherigen Übung verwendet.
Sie können entweder die gespeicherten Daten aus der vorangegangenen Einheit verwenden oder einen neuen Raumausschnitt bestimmen, herunterladen und bearbeiten. Im Prinzip wird jedoch zuerst die Arbeitsumgebung geladen.
# ---- 0 Projekt Setup ----
# Achtung Pakete müssen evtl. manuell installiert werden
library(envimaR)
#--- Schalter für den Download der sentinel daten
get_sen = FALSE
#--- schalter ob digitalisiert werden muss falls er auf FALSE gesetzt ist werden die
# (zuvor erstellten und gesciherten Daten ) im else Teil der Verzeigung eingelesen
digitize = FALSE
## setzen des aktuellen Projektverzeichnisses (erstellt mit envimaR) als root_folder
#root_folder = find_rstudio_root_file()
root_folder = "~/edu/geoinfo/"
# Einlesen des zuvor erstellten Setup-Skripts
source(file.path(root_folder, "src/functions/000_setup.R"))
nclasses=2
Bitte ergänzen Sie (bei auftreteneden Fehlermeldungen) etwaig fehlende oder defekte Pakete im obigen Setup-Skripts.
Auf der Grundlage der verfügbaren Sentinel Daten sollten nun als Erstes geeignete Datensätze für eine Oberflächenklassifikation identifiziert werden. Hierzu kann der vollständige Datensatz auch vom Kursdatenserver heruntergeladen werden (Bitte beachten Sie dass sie im VPN bzw. UniNetz angemeldet sein müssen). Entpacken Sie diese Daten in das Wurzelverzeichnis des Kursprojekts. Das heisst der dort bereits vorhandene Ordner data wird ersetzt/ergänzt.
Bei näherer Betrachtung der RGB Bilder (RGB432B) zeigt sich, das vier Datensätze aufgrund der Bildqualität und geringen Wolkenbedeckung geeignet zu sein scheinen. Es handelt sich um den 19.06. bzw. 24.7. 2019 und den 33.06 bzw. 30.7. 2020. Aufgrund des früheren Vegetationszeitpunktes wurden die Maibilder gewählt, da hier evtl. Aufwuchs auf den Rodungsflächen weniger stark sichtbar sein könnte.
Zunächst müssen diese Daten in einem “Rasterstapel” also einem Mehrkanalbild verfügbar sein. Bereits bei der ersten Beschäftigung mit dem sen2r Paket hatten wir weitere Produkte als sinnvoll erachtet und berechnet. Es sind folgende Indizes:
Diese Daten müssen für den gewünschten Zeitslot und Raumausschnitt besorgt und vorbereitet werden. Aus Gründen der Vereinfachung gibt es für den Kurs eine entsprechende sen2R-Projektdatei:
{
"preprocess": [true],
"s2_levels": ["l1c", "l2a"],
"sel_sensor": ["s2a", "s2b"],
"online": [true],
"server": ["scihub"],
"order_lta": [true],
"downloader": ["builtin"],
"overwrite_safe": [false],
"rm_safe": ["yes"],
"max_cloud_safe": [5],
"step_atmcorr": ["auto"],
"sen2cor_use_dem": [true],
"sen2cor_gipp": [null],
"timewindow": ["2019-06-01", "2020-08-31"],
"timeperiod": ["seasonal"],
"extent": ["{\n \"type\": \"FeatureCollection\",\n \"features\": [\n {\n \"type\": \"Feature\",\n \"properties\": {\n \"X_leaflet_id\": 2080,\n \"layerId\": \"2080\",\n \"edit_id\": \"2080\"\n },\n \"geometry\": {\n \"type\": \"Polygon\",\n \"coordinates\": [\n [\n [10.198975, 51.848575],\n [10.198975, 51.947721],\n [10.413895, 51.947721],\n [10.413895, 51.848575],\n [10.198975, 51.848575]\n ]\n ]\n }\n }\n ]\n}"],
"s2tiles_selected": ["32UNC"],
"s2orbits_selected": [null],
"list_prods": ["BOA", "SCL", "CLD"],
"list_indices": ["ARVI", "EVI", "MSAVI2", "MSI", "NDVI", "SAVI"],
"list_rgb": ["RGB432B", "RGB843B"],
"rgb_ranges": [
[null],
[
[0, 7500],
[0, 2500],
[0, 2500]
]
],
"index_source": ["BOA"],
"mask_type": ["nodata"],
"max_mask": [5],
"mask_smooth": [0],
"mask_buffer": [0],
"clip_on_extent": [true],
"extent_as_mask": [true],
"extent_name": ["harz"],
"reference_path": [null],
"res": [null],
"res_s2": ["10m"],
"unit": ["Meter"],
"proj": [null],
"resampling": ["near"],
"resampling_scl": ["near"],
"outformat": ["GTiff"],
"rgb_outformat": ["GTiff"],
"index_datatype": ["Int16"],
"compression": ["DEFLATE"],
"rgb_compression": ["DEFLATE"],
"overwrite": [true],
"path_l1c": ["~/edu/geoinfo/data/data_lev0"],
"path_l2a": ["~/edu/geoinfo/data/data_lev0"],
"path_tiles": [null],
"path_merged": [null],
"path_out": ["~/edu/geoinfo/data/data_lev1"],
"path_rgb": ["~/edu/geoinfo/data/data_lev1"],
"path_indices": [""],
"path_subdirs": [true],
"thumbnails": [true],
"log": ["~/edu/geoinfo//doc/nw-harz.log"],
"parallel": [true],
"processing_order": ["by_date"],
"pkg_version": ["1.5.0.9000"]
}
#--- Download der Daten
# gui = TRUE ruft die GUI zur Kontrolle auf
if (get_sen){
out_paths_3 <- sen2r(
gui = T,
param_list = "~/edu/geoinfo/data/harz.json",
tmpdir = envrmt$path_tmp,
)
}
#--- Einlesen der Daten aus den Verzeichnissen
# RGB stack der beiden Jahre
pred_stack_2019 = raster::stack(list.files(file.path(envrmt$path_data_lev1,"RGB843B"),pattern = "20190619",full.names = TRUE))
pred_stack_2020 = raster::stack(list.files(file.path(envrmt$path_data_lev1,"RGB843B"),pattern = "20200623",full.names = TRUE))
# Stack-Loop über die Daten
for (pat in c("RGB432B","EVI","MSAVI2","NDVI","SAVI")){
pred_stack_2019 = raster::stack(pred_stack_2019,raster::stack(list.files(file.path(envrmt$path_data_lev1,pat),pattern = "20190619",full.names = TRUE)))
pred_stack_2020 = raster::stack(pred_stack_2020,raster::stack(list.files(file.path(envrmt$path_data_lev1,pat),pattern = "20200623",full.names = TRUE)))
}
# Zuweisen von leserlichen Namen auf die Datenebenen
names(pred_stack_2019) = c("nir","red","green","red","green","blue","EVI","MSAVI2","NDVI","SAVI")
names(pred_stack_2020) = c("nir","red","green","red","green","blue","EVI","MSAVI2","NDVI","SAVI")
saveRDS(pred_stack_2019,paste0(envrmt$path_data,"pred_stack_2019.rds"))
saveRDS(pred_stack_2020,paste0(envrmt$path_data,"pred_stack_2020.rds"))
pred_stack_2019 = readRDS(paste0(envrmt$path_data,"pred_stack_2019.rds"))
pred_stack_2020 = readRDS(paste0(envrmt$path_data,"pred_stack_2020.rds"))
names(pred_stack_2019) = c("nir","red_1","green_1","red_2","green_2","blue","EVI","MSAVI2","NDVI","SAVI")
names(pred_stack_2020) = c("nir","red_1","green_1","red_2","green_2","blue","EVI","MSAVI2","NDVI","SAVI")
# visuelle Überprüfung der stacks
plot(pred_stack_2019)

plot(pred_stack_2020)

# Einlesen der Corine Daten
# Für den Download ist ein Konto notwendig https://land.copernicus.eu/pan-european/corine-land-cover
# Daher die Daten manuell herunterladen und in das Verzeichnis kopieren und entpacken
# Dann auskommentierten Tail ausführen
# corine_eu = raster(file.path(envrmt$path_data_lev0,"u2018_clc2018_v2020_20u1_raster100m/DATA/U2018_CLC2018_V2020_20u1.tif"))
# tmp = projectRaster(pred_stack_2019[[1]],crs = crs(corine_eu))
# corine_crop = raster::crop(corine_eu,tmp)
# corine_utm = projectRaster(corine_crop,crs = crs(pred_stack_2019))
# corine = resample(corine_utm,pred_stack_2019[[1]])
# raster::writeRaster(corine,file.path(envrmt$path_data_lev0,"/corine.tif"),overwrite=TRUE)
# ansonsten den Beipieldatensatz laden
corine = raster::raster(file.path(envrmt$path_data_lev0,"corine.tif"))
# Erstellen einer Wald-Maske
# Agro-forestry areas code=22, Broad-leaved forest code=23,
# Coniferous forest code=24, Mixed forest code=25
mask = reclassify(corine,c(-100,22,0,22,26,1,26,500,0))
mapview(mask)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 1693870
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 1693870 '
Die wohl bekannteste unüberwachte Klassifizierungstechnik ist das K-means-Clustering, das auch als der “einfachster Algorithmus des maschinellen Lernens” bezeichnet wird. Er wird häufig angewendet um einen ersten Überblick zu erhaltne ob die Rasterdaten im Merkmalsraum ausreichend trennbar sind.
In unserem Beispiel (auf 2 Klassen angewendet und mit unsuperClass Funktion aus dem RStoolbox Paket ausgeführt) sieht das wie folgt aus:
## k-means über RStoolbox
# Modell
prediction_kmeans_2019 = unsuperClass(pred_stack_2019, nClasses = 2,norm = TRUE, algorithm = "MacQueen")
# Klassifikation
mapview(prediction_kmeans_2019$map, col = c('darkgreen', 'burlywood', 'green'))
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 1693870
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 1693870 '
prediction_kmeans_2020 = unsuperClass(pred_stack_2020, nClasses = 2,norm = TRUE, algorithm = "MacQueen")
mapview(prediction_kmeans_2020$map, col = c('darkgreen', 'burlywood', 'green'))
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 1693870
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 1693870 '
Trainingsdaten werden vielfach durch abdigitalisieren manuell erfasst. Dies kann recht komfortabel in RStudio durchgeführt werden, falls nur schnell und effektiv einige wenige Trainingsflächen zu digitalisieren sind.
Für größere Aufgaben ist es sinnvoll, auf den hohen Komfort z.b. der QGIS-Arbeitsumgebung zurückzugreifen. Für diese Übung verwenden wir mapedit, ein kleines, aber feines Paket, das die Digitalisieren und Editieren von Vektordaten im Rstudio- oder externen -browser ermöglicht. In Kombination mit mapview können auch beliebige Farbkomposita als Grundlage für die Digitalisierung verwendet werden.
Wir nehmen an, dass wir zwei Arten von Landbedeckung klassifizieren wollen: clearcut (Abholzungen) und other (Anderes). Mit mapedit `muss jede Klasse einzeln digitalisiert werden. Sobald die Trainingsgebiete als Vektordaten verfügbar sind können die Merkmale des jeweiligen Rasterstacks entsprechend der digitalisierten Klassen in eine Tabelle extrahiert und auf etwaige Fehlwerte bereiningt werden.
Falls dieser Teil bereits absolviert wurde kann die logische Variable (am Anfang des Scripts definiert) digitize auf FALSE gesetzt werden und es wird dann der else Teil der Verzweigung durchlaufen - also nur noch die existierenden Daten eingelesen.
Exkurs: Trainingsgebiete mit mapedit erfassen
m1 = tm_shape(pred_stack_2019) + tm_rgb(r=1, g=2, b=3) +
tm_layout(legend.outside.position = "right",
legend.outside = T,
panel.label.height=0.6,
panel.label.size=0.6,
panel.labels = c("r=1, g=2, b=3")) +
tm_grid()
m2 = tm_shape(pred_stack_2019) + tm_rgb(r=4, g=5, b=6) +
tm_layout(legend.outside.position = "right",
legend.outside = T,
panel.label.height=0.6,
panel.label.size=0.6,
panel.labels = c("r=4, g=5, b=6")) +
tm_grid()
tmap::tmap_arrange(m1,m2)
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## Warning: Raster values found that are outside the range [0, 255]
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## Warning: Raster values found that are outside the range [0, 255]
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## Warning: Raster values found that are outside the range [0, 255]
## stars object downsampled to 1151 by 868 cells. See tm_shape manual (argument raster.downsample)
## Warning: Raster values found that are outside the range [0, 255]

Verwenden Sie die Ebenensteuerung, um die Ebenen umzuschalten. Bei Echtfarbkompositen werden die sichtbaren Spektralkanäle Rot (B04), Grün (B03) und Blau (B02) den entsprechenden roten, grünen bzw. blauen Farbkanälen zugeordnet, wodurch ein quasi natürliches “farbiges” Bild der Oberfläche entsteht, wie es ein Mensch sehen würde, der auf dem Satelliten sitzt. Falschfarbenbilder werden häufig mit den Spektralkanälen für das nahe Infrarot, Rot und Grün erzeugt. Sie eignen sich hervorragend für die Einschätzung der Vegetation, da Pflanzen nahes Infrarot und grünes Licht reflektieren, während sie rotes Licht absorbieren (Red Edge Effect). Ein dichterer Pflanzenbewuchs ist dunkler rot. Städte und offener Boden sind grau oder hellbraun, und Wasser erscheint blau oder schwarz.
#---- Digitalisierung der Trainingsdaten ----
if (digitize) {
# Für die überwachte Klassifikation benötigen wir Trainingsgebiete. Sie können Sie wie nachfolgend digitalisieren oder alternativ z.B. QGis verwenden
#--- 2019
# Kahlschlag
# Es wird das Falschfarbenkomosit in originaler Auflösung genutzt (maxpixels = 1693870)
# Bitte beachten Sie dass es (1) deutlich länger lädt und (2) Vegetation in Rot dargestellt wird.
# Die Kahlschäge sind jetzt grün
train_area <- mapview::viewRGB(pred_stack_2019, r = 1, g = 2, b = 3, maxpixels = 1693870) %>% mapedit::editMap()
# Hinzufügen der Attribute class (text) und id (integer)
clearcut <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "clearcut", id = 1)
# other: hier gilt es möglichst verteilt übers Bild möglichst alle nicht zu Kahlschlag gehörenden Flächen zu erfassen.
train_area <- mapview::viewRGB(pred_stack_2019, r = 1, g = 2, b = 3) %>% mapedit::editMap()
other <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "other", id = 2)
# rbind kopiert die beiden obigen Vektorobjekte in eine Datei
train_areas_2019 <- rbind(clearcut, other)
# Umprojizieren auf die Raster Datei
train_areas_2019 = sf::st_transform(train_areas_2019,crs = sf::st_crs(pred_stack_2019))
mapview(train_areas_2019, zcol="class")
# Extraktion der Trainingsdaten für die digitalisierten Flächen
tDF_2019 = exactextractr::exact_extract(pred_stack_2019, train_areas_2019, force_df = TRUE,
include_cell = TRUE,include_xy = TRUE,full_colnames = TRUE,include_cols = "class")
# auch hier wieder zusamenkopieren in eine Datei
tDF_2019 = dplyr::bind_rows(tDF_2019)
# Löschen von etwaigen Zeilen die NA (no data) Werte enthalten
tDF_2019 = tDF_2019[complete.cases(tDF_2019) ,]
tDF_2019 = tDF_2019[ ,rowSums(is.na(tDF_2019)) == 0]
# check der extrahierten Daten
summary(tDF_2019)
mapview(train_areas_2019)+pred_stack_2019[[1]]
# Abspeichern als R-internes Datenformat
# ist im Repo hinterlegt und kann desahlb (zeile drunter) eingeladen werden
saveRDS(tDF_2019, paste0(envrmt$path_data,"train_areas_2019.rds"))
# # ---- Das gleiche muss für 2020 wiederholt werden zum digitalisieren und extrahieren bitte ent-kommentieren ----
# Kahlschlag
train_area <- mapview::viewRGB(pred_stack_2020, r = 1, g = 2, b = 3,maxpixels = 1693870) %>% mapedit::editMap()
clearcut <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "clearcut", id = 1)
train_area <- mapview::viewRGB(pred_stack_2020, r = 1, g = 2, b = 3) %>% mapedit::editMap()
other <- train_area$finished$geometry %>% st_sf() %>% mutate(class = "other", id = 2)
train_areas_2020 <- rbind(clearcut, other)
train_areas_2020 = sf::st_transform(train_areas_2020,crs = sf::st_crs(pred_stack_2020))
mapview(train_areas_2020, zcol="class")
tDF_2020 = exactextractr::exact_extract(pred_stack_2020, train_areas_2020, force_df = TRUE,
include_cell = TRUE,include_xy = TRUE,full_colnames = TRUE,include_cols = "class")
tDF_2020 = dplyr::bind_rows(tDF_2020)
tDF_2020 = tDF_2020[ rowSums(is.na(tDF_2020)) == 0,]
saveRDS(tDF_2020, paste0(envrmt$path_data,"train_areas_2020.rds"))
} else {
train_areas_2019 = readRDS(paste0(envrmt$path_data,"train_areas_2019.rds"))
train_areas_2020 = readRDS(paste0(envrmt$path_data,"train_areas_2020.rds"))
}
Als Resultat ligen nun zwei tabellarische Trainingsdatensätze für 2019 und 2020 vor. Beide Datensätze beinhalten alle Rasterinformationen die von den Polygonen für die Klassen “clearcut” und “other” abgedeckt werden.
head(train_areas_2019)
## class nir red_1 green_1 red_2 green_2 blue EVI MSAVI2 NDVI SAVI X
## 1 clearcut 74 64 58 64 58 39 2969 2648 5524 2978 586485
## 2 clearcut 67 55 49 55 49 32 2778 2498 5684 2854 586495
## 3 clearcut 75 96 73 96 73 49 2191 2018 3975 2300 586475
## 4 clearcut 75 109 78 109 78 57 1962 1792 3464 2049 586485
## 5 clearcut 74 97 73 97 73 53 2161 1954 3876 2234 586495
## 6 clearcut 62 72 56 72 56 38 2100 1897 4422 2221 586505
## Y cell coverage_fraction
## 1 5752735 514566 0.27757052
## 2 5752735 514567 0.08858948
## 3 5752725 516064 0.32311472
## 4 5752725 516065 0.99449599
## 5 5752725 516066 0.98001540
## 6 5752725 516067 0.68315113
head(train_areas_2020)
## class nir red_1 green_1 red_2 green_2 blue EVI MSAVI2 NDVI SAVI X
## 1 clearcut 86 95 72 95 72 53 2793 2542 4606 2820 586645
## 2 clearcut 83 102 74 102 74 52 2475 2297 4190 2569 586655
## 3 clearcut 57 83 54 83 54 39 1576 1453 3473 1732 586665
## 4 clearcut 85 99 72 99 72 50 2603 2427 4390 2701 586645
## 5 clearcut 85 100 71 100 71 51 2615 2435 4389 2708 586655
## 6 clearcut 75 99 68 99 68 48 2154 2007 3915 2285 586665
## Y cell coverage_fraction
## 1 5753385 417147 0.249708712
## 2 5753385 417148 0.474895865
## 3 5753385 417149 0.006561131
## 4 5753375 418646 0.750996768
## 5 5753375 418647 1.000000000
## 6 5753375 418648 0.717114925
Klassifikatoren (z.B der Maximum-Likelihood Klassikator) oder Algorithmen des maschinellen Lernens (wie z. B. Random Forest) ermitteln auf Basis der Trainingsdaten Beschreibungsmodelle die z.b. statistische Signaturen, Klassifikationsbäume oder andere Funktionen darstellen. Im Rahmen der Güte der Trainingsdaten sind solche Modelle geeignet und repräsentativ um Vorhersagen für Räume zu treffen falls die Prädiktoren aus dem Modell flächendeckend vorhanden sind.
Wir wollen nun die räumlichen Merkmale Kahlschlag/kein Wald exemplarisch mit einer Maximum Likelihood Klassifikation und mit Random Forest vorhersagen und Standardmethoden der Zufallsvalidierung und Modellgüteeinschätzung anwenden.
Ziel ist es Kahlschläge von allen anderen Pixeln zu trennen und die Unterschide von 2019/2020 zu quantifizieren.
Exkurs: Random Forest
Random Forests können sowohl für Regressions- als auch für Klassifizierungsaufgaben verwendet werden, wobei letztere besonders in der Umwelt-Fernerkundung relevant sind. Wie bei jeder maschinellen Lernmethode lernt auch das Random-Forest-Modell, Muster und Strukturen in den Daten selbst zu erkennen.
!
Abbildung: Vereinfachte Darstellung der Klassifizierung von Daten durch Random Forest während des Trainings. Venkata Jagannath [CC BY-SA 4.0] via wikipedia.org
Ein Random-Forest-Algorithmus lernt über die Daten, indem er viele Entscheidungsbäume erstellt - daher auch der Name “Forest.” Für Klassifizierungsaufgaben nimmt der Algorithmus eine passende Instanz aus eines Baumes aus dem Trainingsdatensatz und weist dem Pixel dies Klasse zu. Dies wird mit allen verfügbaren Entschidungsbäuumen so gemacht und letzlich wird nach dem Winner takes it all Ansatz die Instanz der Klasse zugeordnet, die das Ergebnis der meisten Bäume ist.
Da der Random-Forest-Algorithmus Trainingsdaten benötigt, handelt es sich um eine überwachte Lernmethode. Das bedeutet, dass wir als Nutzer dem Algorithmus Daten zur Verfügung stellen müssen, die im Wissen über die vorherzusagenden Klassen vermitteln. Diese Daten werden dann in Trainings- und Testdaten aufgeteilt.
## ---- Überwachte mit Klassifikation Random Forest ----
train_areas_2019 = readRDS(paste0(envrmt$path_data,"train_areas_2019.rds"))
train_areas_2020 = readRDS(paste0(envrmt$path_data,"train_areas_2020.rds"))
## Hier wird der Random Forest über das Utility Paket caret aufgerufen
# Setzen eines "seed" ermöglicht reproduzierbaren Zufall
set.seed(123)
# Zufälliges Ziehen von 25% der Daten (Training/Test)
idx_2019 = createDataPartition(train_areas_2019$class,list = FALSE,p = 0.25)
trainDat_2019 = train_areas_2019[idx_2019,]
testDat_2019 = train_areas_2019[-idx_2019,]
idx_2020 = createDataPartition(train_areas_2020$class,list = FALSE,p = 0.25)
trainDat_2020 = train_areas_2020[idx_2020,]
testDat_2020 = train_areas_2020[-idx_2020,]
# Response-Variable (=Spalte "class") muss den Datentyp "factor" haben
trainDat_2019$class <- as.factor(trainDat_2019$class)
testDat_2019$class <- as.factor(testDat_2019$class)
trainDat_2020$class <- as.factor(trainDat_2020$class)
testDat_2020$class <- as.factor(testDat_2020$class)
# Einstellungen Modelltraining: cross-validation, 10 Wiederholungen
ctrlh = trainControl(method = "cv",
number = 10,
savePredictions = TRUE)
#--- random forest model training
cv_model_2019 = train(trainDat_2019[,2:11], # in den Spalten 2 bis 20 stehen die Trainingsdaten (Prediktoren genannt)
trainDat_2019[,1], # in der Spalte 1 steht die zu klassizierende Variable (Response genannt)
method = "rf", # Methode hier rf für random forest
metric = "Kappa", # Qualitäts/Performanzmaß KAppa
trControl = ctrlh, # obig erzeugte Trainingssteuerung soll eingelsen werden
importance = TRUE) # Die Bedeung der Variablen wird mit abgespeichert
cv_model_2019
## Random Forest
##
## 5996 samples
## 10 predictor
## 2 classes: 'clearcut', 'other'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 5397, 5396, 5397, 5396, 5396, 5396, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.9879911 0.8578546
## 6 0.9884903 0.8658856
## 10 0.9879903 0.8600513
##
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 6.
Die Testdaten werden nun für die unabhängige Qualitätsprüfung des Modells verwendet. Eine Wahrheits- oder Konfusionsmatrix zeigt an, wie präzise das Modell die korrekten Klassen vorhersagt. Die Hauptdiagonale der Matrix zeigt die Fälle an, in denen das Modell korrekt ist. In diesem Falle mit nur 2 Klassen gilt der Sonderfall der Beurteilung eines binären Klassifikators. Speziell für die hier verwendete Funktion finden sich ausführliche Erläuertungen in der caret Hilfe.
Die wesentlichen Ausagen zur Modellgüte sind : * ‘Positive’ Class = clearcut: wird mit der Sensitivität (true positive rate) erfasst, die die Wahrscheinlichkeit angibt, mit der ein positives Objekt korrekt als positiv klassifiziert wird. * ‘Negative Class’ = other: wird mit der Spezifität (true negative rate) erfasst und gibt die Wahrscheinlichkeit an, mit der ein negatives Objekt korrekt als negativ klassifiziert wird * Positive and negative predictive values geben für clearcut bzw. other die reale Perfomance an. Sie sind um die reale Häufikeitsverteilung korrigiert und ein Schätzmass für die Präzision bzw. Pereformance des Modells bezüglich der jeweiligen Klassen.
Trotz der hohen Werte sehen wir das die Klasse clearcut hier deutlich abfällt also durchaus noch Verbesserungsbedarf für eine guten Klassifikation vorliegt.
Insgesamt kann das Modell jedoch als gut bezeichnet werden.
# ---- Berechnung der Konfusionsmatrix ----
cm_rf <- confusionMatrix(data = predict(cv_model_2019, newdata = testDat_2019), testDat_2019$class)
cm_rf
## Confusion Matrix and Statistics
##
## Reference
## Prediction clearcut other
## clearcut 699 103
## other 120 17063
##
## Accuracy : 0.9876
## 95% CI : (0.9859, 0.9892)
## No Information Rate : 0.9545
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.8559
##
## Mcnemar's Test P-Value : 0.284
##
## Sensitivity : 0.85348
## Specificity : 0.99400
## Pos Pred Value : 0.87157
## Neg Pred Value : 0.99302
## Prevalence : 0.04554
## Detection Rate : 0.03887
## Detection Prevalence : 0.04459
## Balanced Accuracy : 0.92374
##
## 'Positive' Class : clearcut
##
Nun sind wir soweit das überprüfte Modell auf unseren Datensatz anzuwenden. Üblicherweise wird das in der Fernerkundung Klassifikation genannt.
# Klassifikation (auch Vorhersage genannt)
prediction_rf_2019 = raster::predict(pred_stack_2019 ,cv_model_2019 )
Training und Vorhersage für das Jahr 2020
#--- 2020
# Einstellungen wie 2019
cv_model_2020 = train(trainDat_2020[,2:11],
trainDat_2020[,1],
method = "rf",
metric = "Kappa",
trControl = ctrlh,
importance = TRUE)
cv_model_2020
## Random Forest
##
## 7845 samples
## 10 predictor
## 2 classes: 'clearcut', 'other'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 7061, 7060, 7060, 7060, 7061, 7060, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.9847046 0.9218005
## 6 0.9841944 0.9192613
## 10 0.9838122 0.9174582
##
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.
prediction_rf_2020 = raster::predict(pred_stack_2020 ,cv_model_2020)
Bei der Maximum-Likelihood-Klassifizierung wird davon ausgegangen, dass die Verteilung der Daten für jede Klasse und in jedem Kanal normal verteilt sind. Unter dieser Vorraussetzung wird die Wahrscheinlichkeit berechnet, dass ein bestimmtes Pixel zu einer bestimmten Klasse gehört. Da auch die Wahrscheinlichkeiten als Schwellenwert angegeben werden können werden ohne diese Einschränkung alle Pixel ungeachtet wie unwahrscheinlich zugeordnet. Jedes Pixel wird der Klasse zugeordnet, die die höchste Wahrscheinlichkeit aufweist (d. h. die maximale Wahrscheinlichkeit).

# ---- Maximum Likelihood Classification ----
# superClass() Funktion aus dem Paket RSToolbox umwandeln der Tabelle in das
# geforderte SpatialdataPoint Objekt Identifikation der Koordinaten
sp::coordinates(trainDat_2019) = ~X+Y
sp::coordinates(trainDat_2020) = ~X+Y
# Zuweisen der korrekten Projektion (hier aus dem zugehörigen Raster-Stack)
crs(trainDat_2019) = crs(pred_stack_2019)
crs(trainDat_2020) = crs(pred_stack_2020)
# superClass method "mlc" berechnet die Statisik und klassifiziert in einem aufruf
prediction_mlc_2019 <- superClass(pred_stack_2019, trainData = trainDat_2019[,1:11], responseCol = "class",
model = "mlc", tuneLength = 1, trainPartition = 0.5,verbose = FALSE)
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
prediction_mlc_2020 <- superClass(pred_stack_2020, trainData = trainDat_2020[,1:11], responseCol = "class",
model = "mlc", tuneLength = 1, trainPartition = 0.5,verbose = FALSE)
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
## Warning in method$fit(x = x, y = y, wts = wts, param = tuneValue, lev =
## obsLevels, : Covariance matrix of class/classes clearcut, other is singular,
## i.e. holds perfectly correlated variables.
# Ergebnisse der rf Klassifikation
#prediction_mlc_2019$model
#prediction_mlc_2020$model
## ---- Visualisierung mit mapview ----
mapview::viewRGB(mask*pred_stack_2020, r = 4, g =5, b = 6,maxpixels = 1693870)+
mapview(mask*prediction_rf_2019 , alpha.regions = 0.5, maxpixels = 1693870,
col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = TRUE) +
mapview(mask*prediction_rf_2020, alpha.regions = 0.5, maxpixels = 1693870,
col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
mapview(mask*prediction_mlc_2019$map,alpha.regions = 0.5, maxpixels = 1693870,
col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE) +
mapview(mask*prediction_mlc_2020$map,alpha.regions = 0.5, maxpixels = 1693870,
col.regions = mapviewPalette("mapviewRasterColors"),at = seq(0, nclasses, 1), legend = FALSE)
## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same
Betrachten sie sie nachfolgenden Ressourcen als Beispiele dafür, wie aus solchen Quellen (die alle mehr oder weniger technisch ähnlich sind) ein bestimmtes Set von Werkzeugen zur Bearbeitung einer Fragestellung “entsteht”. Dann nach viel Recherche und kritischer Prüfung wird ein aktueller Stand der Technik innerhalb der wissenschaftlcihen Gemeinschaft erreicht der sich zudem natürlich stetig weiterentwickelt.
Schauen sie einmal in die nachstehende Auswahl von Blogs und Hilfestellungen:
Digitalisierung von Trainingsdaten by the Forest Inventory and Remote Sensing department at the University of Goettingen (Germany)
Digitalisierung Turorial in the QGIS 3.16 documentation
Robert J. Hijmans rspatial - supervised classification
Ivan Lizarazo RPubs Tutorial
Sydney Goldstein blog
João Gonçalves supervised classification
Valentin Stefan pixel-based supervised classification
Die beiden ersten Beiträge sind rein technische Tutorials. Die weiteren Beiträge mischen technische Anleitungen mit konzeptionellen oder konkreten fachlichen Fragestellungen. Dennoch sind sie hier nicht als fachwissenschaftliche Referenz gedacht. Sie dienen vielmehr dazu zu zeigen wie technisches und konzeptionelles Verständnis schrittweise erarbeiet werden kann und unterstützen, durch “nachkochen” und anwenden, eigenständiger die Lösung von Fragestellungen anzugehen.
Ich möchte ausdrücklich Valentin Stefan den Autor des Blogbeitrags pixel-based supervised classification zitieren:
“[…] Betrachten Sie diesen Inhalt als einen Blogbeitrag und nichts weiter. Er erhebt nicht den Anspruch, eine erschöpfende Übung oder ein Ersatz für Ihr kritisches Denken zu sein. […]”